
module mo_photoin

  use shr_kind_mod,     only : r8 => shr_kind_r8
  use cam_logfile,      only : iulog
  use cam_abortutils,   only : endrun

  implicit none

  save

  public :: photoin_inti
  public :: photoin
  private

  integer               :: jo2_ndx = 0

  logical, allocatable  :: z_dep(:)
  character(len=32), allocatable :: pht_tag(:)

contains

  subroutine photoin_inti( nlng, lng_indexer )
    !-------------------------------------------------------------
    ! 	... assign use masks
    !-------------------------------------------------------------

    use mo_params,   only : largest
    use mo_setcld,   only : setcld_inti
    use chem_mods,   only : phtcnt, rxt_tag_lst

    implicit none

    !-------------------------------------------------------------
    ! 	... dummy arguments
    !-------------------------------------------------------------

    integer, intent(in) :: nlng
    integer, intent(in) :: lng_indexer(:)

    !-------------------------------------------------------------
    ! 	... local variables
    !-------------------------------------------------------------
    integer :: astat
    integer :: m
    integer :: ndx
    character(len=32) :: jname

    !-------------------------------------------------------------
    ! 	... allocate module arrays
    !-------------------------------------------------------------
    has_photorates : if( nlng > 0 ) then
       allocate( z_dep(nlng), pht_tag(nlng), stat=astat )
       if( astat /= 0 ) then
          write(iulog,*) 'photoin_inti: failed to allocate z_dep; error = ',astat
          call endrun
       end if
       ndx = 0
       do m = 1,phtcnt
          if( lng_indexer(m) > 0 ) then
             if( any( lng_indexer(:m-1) == lng_indexer(m) ) ) then
                cycle
             end if
             ndx = ndx + 1
             pht_tag(ndx) = trim( rxt_tag_lst(m))
          end if
       end do
       if( ndx /= nlng ) then
          write(iulog,*) 'photoin_inti: corrupted lng_indexer'
          call endrun
       end if
       write(iulog,*) ' '
       write(iulog,*) 'photoin_inti: lng_indexer name mapping'
       write(iulog,'(5a)') pht_tag(:)
       write(iulog,*) ' '
       !-------------------------------------------------------------
       ! 	... search for jo2
       !-------------------------------------------------------------
       do m = 1,nlng
          if( pht_tag(m) == 'jo2' ) then
             jo2_ndx = m
             exit
          end if
       end do
       write(iulog,*) ' '
       write(iulog,*) 'photoin_inti: jo2 index = ',jo2_ndx
       write(iulog,*) ' '
       !-------------------------------------------------------------
       ! 	... set altitude dependence logical array
       !-------------------------------------------------------------
       z_dep(:) = .true.
       do m = 1,nlng
          jname = pht_tag(m)
          select case( jname )
          case( 'jno2', 'jno3', 'jho2', 'jhno2', 'jho2no2' )
             z_dep(m) = .false.
          case( 'jc2h5cho', 'jchocho', 'jch3ooh' )
             z_dep(m) = .false.
          end select
       end do

       !-------------------------------------------------------------
       ! 	... intialize cloud layer module
       !-------------------------------------------------------------
       call setcld_inti

    end if has_photorates

  end subroutine photoin_inti

  subroutine photoin( idate, alat, along, &
       ut, esfact, o3top, o2top, albedo, &
       z, tlev, tlay, xlwc, xfrc, &
       airlev, aocs1, aocs2, acbs1, acbs2, &
       asoa, aant, aso4, asal, ads, o3, rh, &
       prate,  zen, nw, dt_xdiag )
    !----------------------------------------------------------
    !     	... interactive photolysis interface routine
    !----------------------------------------------------------
    
    use mo_tuv_inti, only : nlng
    use mo_params,   only : kj, kw
    use mo_wavelen,  only : deltaw, sflx, wc, wl, wu
    use mo_wavelab , only : sj
    use mo_zadj,     only : adj_coeffs
    use mo_setair,   only : setair
    use mo_setozo,   only : setozo
    use mo_pchem,    only : pchem
    use mo_sphers,   only : sphers
    use mo_airmas,   only : airmas
    use mo_setz,     only : setz
    use mo_seto2,    only : set_o2_xsect
    use mo_rtlink,   only : rtlink
    use mo_setcld,   only : setcld   !, mreg
    use mo_setaer,   only : setaer
    use ppgrid,      only : pver, pverp

    implicit none

    !----------------------------------------------------------
    !     	... dummy arguments
    !----------------------------------------------------------
    integer, intent(in)     ::  idate
    integer, intent(in)     ::  nw
    real(r8), intent(in)    ::  alat, along, o3top, o2top
    real(r8), intent(in)    ::  ut, esfact
    real(r8), intent(in)    ::  zen
    real(r8), intent(in)    ::  albedo(kw)
    real(r8), intent(in)    ::  tlay(pver)
    real(r8), intent(in)    ::  xlwc(pverp)           ! cloud water (g/m3)
    real(r8), intent(in)    ::  xfrc(pverp)           ! cloud fraction
    real(r8), intent(in)    ::  tlev(pverp)
    real(r8), intent(in)    ::  airlev(pverp)
    real(r8), intent(in)    ::  z(pverp)
    real(r8), intent(in)    ::  aocs1(pverp)
    real(r8), intent(in)    ::  aocs2(pverp)
    real(r8), intent(in)    ::  acbs1(pverp)
    real(r8), intent(in)    ::  acbs2(pverp)
    real(r8), intent(in)    ::  asoa(pverp)
    real(r8), intent(in)    ::  aant(pverp)
    real(r8), intent(in)    ::  aso4(pverp)
    real(r8), intent(in)    ::  asal(pverp,4)
    real(r8), intent(in)    ::  ads(4,pverp)
    real(r8), intent(in)    ::  rh(pverp)
    real(r8), intent(inout) ::  o3(pverp)
    real(r8), intent(out)   ::  prate(pverp,nlng)
    real(r8), intent(out)   ::  dt_xdiag(:)

    !----------------------------------------------------------
    !     	... local variables
    !----------------------------------------------------------
    integer  :: i, j, k, km, wn, n, astat
    real(r8) :: factor, delzint
    real(r8) :: wcen
    real(r8), allocatable :: xs(:,:,:)
    real(r8), allocatable :: adjcoe(:,:)    ! ftuv adjustment factor

    !----------------------------------------------------------
    ! 	... altitude grid
    !----------------------------------------------------------
    real(r8)    :: colinc(pverp)
    real(r8)    :: vcol(pverp)
    real(r8)    :: scol(pverp)
    real(r8)    :: to3(pverp)

    !----------------------------------------------------------
    ! 	... solar zenith angle
    !           slant pathlengths in spherical geometry
    !----------------------------------------------------------
    integer     :: nid(0:pver)
    real(r8)    :: dsdh(0:pver,pver)

    !----------------------------------------------------------
    ! 	... extra terrestrial solar flux and earth-sun distance ^-2
    !----------------------------------------------------------
    real(r8)    :: etf(nw)
    real(r8)    :: delw(nw)
    real(r8)    :: xsec(nw)

    !--------------------------------------------------------------
    ! 	... atmospheric optical parameters:
    !--------------------------------------------------------------
    integer, parameter :: mreg = 16   
    integer :: nreg                  ! regions at each grid
    real(r8)    :: dtrl(pver,nw)
    real(r8)    :: dto3(pver,nw)
    real(r8)    :: dto2(pver,nw)
    real(r8)    :: dtcld(pver,nw)
    real(r8)    :: omcld(pver,nw)
    real(r8)    :: gcld(pver,nw)

    real(r8)    :: dtcbs1(pver,nw)
    real(r8)    :: dtcbs2(pver,nw)
    real(r8)    :: omcbs1(pver,nw)
    real(r8)    :: omcbs2(pver,nw)
    real(r8)    :: gcbs1(pver,nw)
    real(r8)    :: gcbs2(pver,nw)

    real(r8)    :: dtocs1(pver,nw)
    real(r8)    :: dtocs2(pver,nw)
    real(r8)    :: omocs1(pver,nw)
    real(r8)    :: omocs2(pver,nw)
    real(r8)    :: gocs1(pver,nw)
    real(r8)    :: gocs2(pver,nw)

    real(r8)    :: dtant(pver,nw)
    real(r8)    :: omant(pver,nw)
    real(r8)    :: gant(pver,nw)

    real(r8)    :: dtsoa(pver,nw)
    real(r8)    :: dtso4(pver,nw)
    real(r8)    :: omso4(pver,nw)
    real(r8)    :: gso4(pver,nw)

    real(r8)    :: dtsal(pver,nw,4)
    real(r8)    :: omsal(pver,nw,4)
    real(r8)    :: gsal(pver,nw,4)

    real(r8)    :: dtds1(pver,nw)
    real(r8)    :: dtds2(pver,nw)
    real(r8)    :: dtds3(pver,nw)
    real(r8)    :: dtds4(pver,nw)
    real(r8)    :: omds1(pver,nw)
    real(r8)    :: omds2(pver,nw)
    real(r8)    :: omds3(pver,nw)
    real(r8)    :: omds4(pver,nw)
    real(r8)    :: gds1(pver,nw)
    real(r8)    :: gds2(pver,nw)
    real(r8)    :: gds3(pver,nw)
    real(r8)    :: gds4(pver,nw)

    real(r8)    :: optr(pver,mreg)         ! cld opt (z dependent) at each region
    real(r8)    :: fp(mreg)                ! probability at each region

    real(r8)    :: xso2(nw,pverp)

    !--------------------------------------------------------------
    ! 	... spectral irradiance and actinic flux (scalar irradiance):
    !--------------------------------------------------------------
    real(r8)    :: radfld(pverp,nw)
    real(r8)    :: radxx(pverp,nw)

    !-------------------------------------------------------------
    !  	... j-values:
    !-------------------------------------------------------------
    integer :: jn, m
    
    !-------------------------------------------------------------
    ! 	... location and time
    !-------------------------------------------------------------
    integer  :: iyear, imonth, iday
    real(r8) :: dtime, ut0

    !-------------------------------------------------------------
    ! 	... allocate wrking xsection array
    !-------------------------------------------------------------
    allocate( xs(nw,pverp,nlng), adjcoe(pverp,nlng), stat=astat )
    if( astat /= 0 ) then
       write(iulog,*) 'photoin: failed to allocate xs, adjcoe; error = ',astat
       call endrun
    end if

    etf(1:nw) = sflx(1:nw) * esfact   ! earth-sun distance effect
    !-------------------------------------------------------
    !  	... air profile and rayleigh optical depths (inter-face)
    !-------------------------------------------------------
    call setair( z, nw, wc, airlev, dtrl, colinc, o2top )

    !-------------------------------------------------------------
    ! 	... ozone optical depths (must give temperature) (inter-face)
    !-------------------------------------------------------------
    call setozo( z, nw, wl, tlay, dto3, to3, o3, airlev, o3top )

    !-------------------------------------------------------------
    ! 	... cloud optical depths
    !-------------------------------------------------------------
    call setcld( z, xlwc, xfrc, nreg, fp, optr )

    !-------------------------------------------------------------
    ! 	... aerosol optical depths
    !-------------------------------------------------------------
    call setaer( z, airlev, rh, aocs1, aocs2, &
         acbs1, acbs2,&
         aant, aso4, asal, ads, asoa, &
         dtcbs1, dtcbs2, omcbs1, omcbs2, gcbs1, gcbs2, &
         dtocs1, dtocs2, omocs1, omocs2, gocs1, gocs2, &
         dtant, omant, gant, &
         dtso4, omso4, gso4, &
         dtsal, omsal, gsal, &
         dtds1, dtds2, dtds3, dtds4, &
         omds1, omds2, omds3, omds4, &
         gds1, gds2, gds3, gds4, dtsoa, nw )
    dt_xdiag(1) = sum( dtcbs1(:,16) + dtcbs2(:,16) )
    dt_xdiag(2) = sum( dtocs1(:,16) + dtocs2(:,16) )
    dt_xdiag(3) = sum( dtso4(:,16) )
    dt_xdiag(4) = sum( dtant(:,16) )
    dt_xdiag(5) = sum( dtsal(:,16,1) + dtsal(:,16,2) + dtsal(:,16,3) + dtsal(:,16,4) )
    dt_xdiag(6) = sum( dtds1(:,16) + dtds2(:,16) + dtds3(:,16) + dtds4(:,16) )
    dt_xdiag(7) = sum( dtsoa(:,16) )
    dt_xdiag(8) = sum( dt_xdiag(1:6) )

    !------------------------------------------------------------
    ! 	... photo-chemical and photo-biological weigting functions. 
    !           for pchem, need to know temperature and pressure profiles.
    !           output:
    !           from pchem:  sj(kj,kz,kw) - for each reaction
    !-------------------------------------------------------------
    xs(:,:,1:nlng) = sj(:,:,1:nlng)
    call pchem( nw, wl, wc, tlev, &
         airlev, nlng, pht_tag, xs )

    !-------------------------------------------------------------
    ! 	... slant path lengths for spherical geometry
    !-------------------------------------------------------------
    call sphers( z, zen, dsdh, nid )
    call airmas( z, zen, dsdh, nid, colinc, vcol, scol )

    !---------------------------------------------------------------
    !    	... modification of coefficent of j-vales function of to3 and zenith
    !---------------------------------------------------------------
    call setz( to3, tlev, adj_coeffs, zen, adjcoe, pht_tag )

    !------------------------------------------------------------------
    ! 	... effective o2 optical depth (sr bands, must know zenith angle!)
    !           assign o2 cross section to sj(1,*,*)
    !------------------------------------------------------------------
    call set_o2_xsect( z, nw, wl, colinc, vcol, scol, dto2, xso2 )
    if( jo2_ndx > 0 ) then
       xs(:,:,jo2_ndx) = xso2(:,:)
    end if

    delw(:nw) = deltaw(:nw) * etf(:nw)

    !---------------------------------------------------
    !  	... monochromatic radiative transfer:
    !           outputs are  fdir, fdn, fup
    !---------------------------------------------------

    ! set for cloud only

    do wn = 1,nw
       radfld(:,wn) = 0._r8 
       omcld(:,wn)  = .9999_r8
       gcld (:,wn)  = .85_r8
    end do

    Cld_reg_loop : do n = 1,nreg
       factor = fp(n)
       do wn = 1,nw
          dtcld(:,wn) = optr(:,n)
       end do

#ifdef NO_AEROSOL
       dtcbs1(:,:)  = 0._r8
       dtcbs2(:,:)  = 0._r8
       dtocs1(:,:)  = 0._r8
       dtocs2(:,:)  = 0._r8
       dtant(:,:)   = 0._r8
       dtso4(:,:)   = 0._r8
       dtsal(:,:,:) = 0._r8
       dtds1(:,:)   = 0._r8
       dtds2(:,:)   = 0._r8
       dtds3(:,:)   = 0._r8
       dtds4(:,:)   = 0._r8

       omcbs1(:,:)  = 0._r8
       omcbs2(:,:)  = 0._r8
       omocs1(:,:)  = 0._r8
       omocs2(:,:)  = 0._r8
       omant(:,:)   = 0._r8
       omso4(:,:)   = 0._r8
       omsal(:,:,:) = 0._r8
       omds1(:,:)   = 0._r8
       omds2(:,:)   = 0._r8
       omds3(:,:)   = 0._r8
       omds4(:,:)   = 0._r8

       gcbs1(:,:)  = 0._r8
       gcbs2(:,:)  = 0._r8
       gocs1(:,:)  = 0._r8
       gocs2(:,:)  = 0._r8
       gant(:,:)   = 0._r8
       gso4(:,:)   = 0._r8
       gsal(:,:,:) = 0._r8
       gds1(:,:)   = 0._r8
       gds2(:,:)   = 0._r8
       gds3(:,:)   = 0._r8
       gds4(:,:)   = 0._r8
#endif
       call rtlink( z, nw, albedo, zen, dsdh, &
            nid, dtrl, dto3, dto2, & 
            dtcld, omcld, gcld, &
            dtcbs1, omcbs1, gcbs1, &
            dtcbs2, omcbs2, gcbs2, &
            dtocs1, omocs1, gocs1, &
            dtocs2, omocs2, gocs2, &
            dtant, omant, gant, &
            dtso4, omso4, gso4, &
            dtsal, omsal, gsal, &
            dtds1, omds1, gds1, &
            dtds2, omds2, gds2, &
            dtds3, omds3, gds3, &
            dtds4, omds4, gds4, radxx )
       do wn = 1,nw
          radfld(:,wn) = radfld(:,wn) + radxx(:,wn)*factor
       end do
    end do Cld_reg_loop

    !----------------------------------------------------------
    !     	... interplation at the top level
    !----------------------------------------------------------
    delzint = (z(pver-1) - z(pver-2))/(z(pver) - z(pver-1))
    do wn = 1,nw
       radfld(1,wn) = radfld(2,wn) + (radfld(2,wn) - radfld(3,wn))*delzint
       radfld(1,wn) = max( radfld(1,wn),radfld(2,wn) )
    end do

    !----------------------------------------------------------
    !   	... j-val calculation
    !           spherical irradiance (actinic flux)
    !           as a function of altitude
    !           convert to quanta s-1 nm-1 cm-2
    !           (1.e-4 * (wc*1e-9) / (hc = 6.62e-34 * 2.998e8))
    !----------------------------------------------------------
    rate_loop : do m = 1,nlng
       if( .not. z_dep(m) ) then
          xsec(:nw)     = xs(:nw,1,m) * delw(:nw)
          prate(:pverp,m) = matmul( radfld, xsec )
       else
          do k = 1,pverp
             km = pverp - k + 1
             xsec(:nw) = xs(:nw,km,m) * delw(:nw)
             prate(k,m) = dot_product( radfld(k,:nw), xsec(:nw) )
          end do
       end if
       prate(1:pverp,m) = prate(1:pverp,m) * adjcoe(pverp:1:-1,m)  
    end do rate_loop

    deallocate( xs, adjcoe )

  end subroutine photoin

end module mo_photoin
